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We explore the nature of the transition to the Fulde-Ferrell-Larkin- Ovchinnikov superfluid phases 
in the low temperature range in two dimensions, for the simplest isotropic BCS model. This is done 
by applying the Larkin-Ovchinnikov approach to this second order transition. We show that there is 
a succession of transitions toward ever more complex order parameters when the temperature goes 
to zero. This gives rise to a cascade with, in principle, an infinite number of transitions. Except for 
one case, the order parameter at the transition is a real superposition of cosines with equal weights. 
The directions of these wavevectors are equally spaced angularly, with a spacing which goes to zero 
when the temperature goes to zero. This singular behaviour in this T — limit is deeply linked to 
the two-dimensional nature of the problem. 

PACS numbers : 74.20.Fg, 74.60.Ec 



I. INTRODUCTION 

The possible existence of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superfluid phases [1,2] has been pointed out 
in the early sixties and it has given rise to much work ever since that time. In addition to their intrinsic fundamental 
interest [3], these phases are quite relevant experimentally since they are expected to arise in superconductors with 
very high critical fields, which are naturally very actively searched for. On several occasions these phases have been 
claimed to be observed experimentally, but to date these hopes have not been firmly substantiated. Very recently 
anomalies in the heavy fermion compound CeCoIns have been attributed to FFLO phases [4]. The case of two- 
dimensional (2D) systems is of particular interest [5] since they are experimentally quite relevant. Indeed a major 
strategy to observe these transitions is to eliminate orbital currents, which are responsible for the low critical fields 
in standard superconductors. This can be achieved in quasi two- dimensional systems, made of widely separated 
conducting planes, such as organic compounds or high T c cuprate superconductors. In this case hopping between 
planes is very severely restricted. Hence the orbital currents perpendicular to the planes are very weak when a strong 
magnetic field is applied parallel to the planes, and there is essentially no orbital pair breaking effect which opens the 
path to FFLO phases at much higher fields. Indeed experimental results in organic compounds have been claimed 
quite recently [6,7] to be compatible with the existence of FFLO phases. Naturally when the magnetic field is not 
exactly parallel to the planes, one finds in addition vortex-like structures and the physical situation gets even more 
complex [8]. 

The FFLO transition in 2D systems is believed to be second order and in particular Burkhardt and Rainer [9] 
have studied in details the transition to a planar phase, where the order parameter A(r) is a simple cos(q.r) at the 
transition. This phase has been found by Larkin and Ovchinnikov [2] to be the best one in 3D at T = for a second 
order phase transition. And in 3D it is also found to be the preferred one in the vicinity of the tricritical point 
and below [10-12], although in this case the transition turns out to be first order (except at very low temperature). 
However it is not clear that this is always the case since, as first explored by Larkin and Ovchinnikov, this order 
parameter is in competition with any superposition of plane waves, provided that their wavevectors have all the same 
modulus. Indeed we have shown very recently [13] that, at low temperature, the transition is rather a first order 
one, toward an order parameter with a more complex structure. For example at T = it is very near the linear 
combination of three cosines oscillating in orthogonal directions. 

In this paper we explore the low temperature range in 2D and show that the second order transition is indeed 
toward rather more complex order parameters. A short report of our results has already been published [14]. A first 
step in this direction is found in the recent work of Shimahara [15] who found a transition toward a superposition of 
three cosines. Here we show that, when the temperature is lowered toward T = 0, one obtains a cascade of transitions 
toward order parameters with an ever increasing number of plane waves. The T = limit is singular in this respect. 
This is actually clear from the beginning. Indeed if one looks at the second order term in the expansion of the free 
energy in powers of the order parameter, which gives the location of the FFLO transition, one finds it to be a singular 
function of the plane wave wavevector. This is recalled in the next section. Then we calculate the fourth order term in 
the free energy expansion and show that the phases which are selected by this term display the cascade of transitions 
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mcntionned above. 



II. THE FREE ENERGY EXPANSION : SECOND ORDER TERM 



The general expression for the free energy difference O = fl s — fl n between the superconducting and the normal 
state can be obtained in a number of ways, starting for example [9,16,17] from Eilcnbcrger's expression in terms of 
the quasiclassical Green's function or from the gap equation [2] and Gorkov's equations. When the result is expanded 
up to fourth order term in powers of the Fourier components A q of the order parameter A(r): 

A(r) =^A qi exp(z qj .r) (1) 

one obtains: 

§-=Yl ° 2 ^ T)\ A q | 2 + °4(qi, Q2, q 3 , q 4 , P, T)A qi A* 2 A qs A* 4 (2) 
q q; 

where we have momentum conservation q! + q 3 = q 2 + q4 in the fourth order term and N is the single spin density 
of states at the Fermi surface. The explicit expression of ^2(9, P, T) in terms of the standard BCS interaction V and 
of the free fermions propagator is: 

7Vo0 2 (<z,M,T) = i-r^G(k)G(k + q) (3) 

where G(k) = (iu> n — £k) _1 and G(k) = (— iu> n — Ck) 1 and tD„ = u> n — ip, with p = (/if — /i|)/2 being half the 
chemical potential difference between the two fermionic populations forming pairs, £k the kinetic energy measured 
from the Fermi surface for p = and ui n — irT(2n + 1) the Matsubara frequency. Performing the £k integration and 
the 2D angular average over k gives: 

n 2 (g, }i, T) = -i- + 2^T Im £ i (4) 

where we have introduced the dimensionless wavevector q = qvp/2p,. In Eq.(4) the summation has to be cut-off at a 
frequency w c in the standard BCS way. It is more convenient to rewrite ^2, by introducing physical quantities related 
to the q — case, as: 

Q 2 (q,fl,T) = ao(fl,T)+I(q,fl,T) (5) 

with: 

00 

I(q, H,T) = 2nT Im £ - — (6) 

We have introduced: 

«o(M,T) = -l--2^rRe£^- (7) 

which is zero on the spinodal transition line (the line in the p,, T plane where the normal state becomes absolutely 
unstable against a transition toward a space independent order parameter) and is positive above it. At the FFLO 
transition we are looking at, we have ^2(9, p, T) = 0. The actual transition corresponds to the largest possible p at 
fixed T. From Eq.(7) this corresponds to have the largest a (p, T). Hence from Eq.(5) we want to minimize I(q, p, T) 
with respect to q. At low temperature it is more convenient to express I as: 

1 f°° uj 1 1 

/(», ft T) . -5*/^ .a„ M -) [ 7 __ - _] (S) 
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where the integration contour runs actually infinitcsimally above the real tu axis. 
At T — the integration is easily performed to give: 



I(q, p, T) = Re ln(l + yj\ - q 2 ) - In 2 (9) 

The minimum is reached for q = 1 , in agreement with Shimahara [5] and Burkhardt and Rainer [9] , and I = — In 2 at 
this minimum. In this case, from Eq.(7), a (p, 0) = ln(2/x/A ) where A = 2to c exp(— 1/N V) is the zero temperature 
BCS phase gap (this corresponds to the value p = Aq/2 for the spinodal transition). This leads to p = Ao for 
the location of the FFLO transition, again in agreement with previous work. It is worth to note that, as already 
mentionned in the introduction, the location q — 1 of the minimum corresponds to a singular point for I since we 
have explicitely / = ln(<j/2) for q > 1. While / itself is continuous, its derivative is discontinuous for q = 1. 

For T / there is no singular behaviour and we find the value of q giving the minimum I by writing that its 
derivative with respect to q is zero. Integrating the result by parts leads to the condition: 

Rc f°° dy —L- - 1 + 2ty = = 2 (10) 
7-00 cosh 2 y ^/(l + 2ty) 2 -q 2 

where we have taken the new variable y = uj/2T and defined the reduced temperature t = T/p. Only the ranges 
y > (q— l)/2t = a and y < —(q + l)/2t contribute to the real part of the integral. Since at low T we have q s=s 1, this 
last range will only give an exponentially small contribution because of the factor cosh -2 y. Since this same factor 
makes \y\ to be at most of order 3 , we can make at low T to leading order l + 2ty ss 1 and (l + 2fy) 2 — q 2 w 2(1 — q+2ty). 
It is then seen that we must have (q — l)/2t » 1, because (q — l)/2t w 1 makes the left hand side of Eq.(10) much 
larger than unity at low T. This implies y » 1 and 2 cosh y w expy. The integral is then easily evaluated and Eq.(10) 
gives finally to leading order: 

9-l = £l*£ (11) 

H 2 2t v ' 

Note that this result is in disagrccmccnt with the analysis given [18] by L. N. Bulaevskii. The reason for this 
discrepancy is discussed in details in Appendix A. In particular we rederive in this appendix our Eq. (10) from the 
starting equation of Ref. [18]. 

A more complete low temperature expansion can be fairly easily extracted from Eq.(10). As previously seen, the 
range y > a in the integration is sufficient since the other integration range gives an exponentially small term. With 
the change of variable y = u 2 + a, Eq.(10) leads to: 



r + co 

L ' 



du 2 \ +2lu ] „ = 2V? (12) 
o cosh (u 2 + a)v 1 + tu 2 
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where we have defined t = t/q. Neglecting terms of order P Eq. (12) can be written as: 

+°° du 3_ f + °° u 2 du = 

+ o* / = 2 ^ (I 3 ) 

cosh (u 2 + a) 2 J cosh (u 2 + a) 

From Eq.(ll) the leading order for a is ao = (1/4) ln(7r/2t) ^> 1. In the second integral in Eq.(13) we can replace a by 
a since exp(— 2a ) = y '2i ~/n <C 1. Setting a = a + 5a, we obtain an expansion of 5a in powers of V-l 2 by expanding 
the first integral up to second order in powers of 5a: 

-^-t+ ^lp' 2 + 5a (-At 1 ' 2 + ^=i] + AiSafi 1 ' 2 + ty 2 = (14) 

valid up to the order P' 2 . Solving this equation order by order, we find the following expansion for a: 

This expression leads to a marked improvement when it is compared to the straight numerical evaluation of the 
Matsubara sums in Eq.(6). This is seen in Fig.l where we have plotted the optimal a = (q — l)/2t from straight 
numerical calculation, as well as its leading low temperature approximation and the one resulting from the expansion 
Eq.(15). We give also for convenience, in the lower panel, the dependence of our reduced temperature t as a function 
of the ratio T/T c0 , where T c0 is the standard BCS critical temperature, found for p = 0. In the low temperature limit 
T — > 0, we have merely q — 1 and p = A , hence t = T/A and the dotted straight line in the lower panel gives this 
limiting behaviour T/T c0 = tA /T c0 . 
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FIG. 1. a = (q — l)/2t from the leading low temperature approximation Eq.(ll) (dashed-dotted line), the low temperature 
expansion Eq. (15) (dotted line) and the full numerical calculation (full line), as a function of i = t/q. The lower panel gives 
the relation between i and T/T c q (full line), the dotted line is the low temperature limiting behaviour T/T c o = iAo/TeO = 1.76?. 

III. FOURTH ORDER TERM 
A. Leading behaviour 

The second order term in the free energy gives the same transition line for all the combinations of plane waves. 
The selection among the various possible order parameters will be made by the fourth order term in the free energy 
expansion, as pointed out by Larkin and Ovchinnikov. The selected state will correspond to the lowest fourth order 
term. So we turn now to this fourth order term. Its general expression is given by Eq.(2) with [2] : 

7Von4(qi,q 2 ,q 3 ,q4,/2,T)=T^G(k)G(k + q 1 )G(k + q 1 -q 4 )G(k + q 2 ) (16) 

71, fc 

where we have used already momentum conservation qi + q3 = q2 + q4- Since all the wavevectors q^ have the same 
modulus given by Eq.(ll), this momentum conservation implies in our 2D situation either [2] qi = q2 together with 
Q3 = Q4, or equivalently qi = q4 together with q3 = q 2 . Moreover one has the additional possibility qi = — q3 
together with q 2 = — q4- The two first possibilities lead to a same coefficient called t/(a qi ,q 3 ) by LO where a qi , q3 
is the angle between qi and q3. Similarly the last possibility leads to a coefficient </(a qi . q2 ). Explicitely the fourth 
order term in the free energy (last term of the r.h.s. of Eq.(2)) becomes [2] : 
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±£(2 - ^, q3 )|A q J 2 |A q3 | 2 J(a q „ qj ) + (1 - 6^ «y qi ,_^)A q4 A_ qi A^A* q /K ( , q J (17) 

with, after a change of variable in the summation: 

^oJ(a qi>qa ) = T^G(k + qi )G(k + q 2 )G 2 (k) (18) 

and 

No J(a qi , q2 ) = t]T G(k)G(k + qi + q 2 )G(k + qi )G(k + q 2 ) (19) 
When the £k integration is performed, taking into account £k+q ~ £k + k.q/m since q <C ftp, one finds: 



CO ^ 

^(a qi .q 2 ) - 4ttT Im V < — — + 1 <-► 2 > (20) 

r ^ (2i^„ + k.qi/m) 2 (2iw„ + k.q 2 /m) 

where the bracket means the angular average, and similarly: 



J Kl , q2 ) = -16^Img < [(2i _ n)2 _ (kqi/m)2 ^ (2i - ri)2 _ (k q2/m)2] > (2D 

on which it is clear that J {it — a) = J (a). 

Let us first consider J (a). As in the preceding section when going from Eq.(6) to Eq.(8), it is more convenient to 
transform the sum over Matsubara frequencies into an integral on the real frequency axis. When one furthermore 
performs a by parts integration over the frequency cj, one gets from Eq.(20): 

/oo r 2n 2 

-oo dUJ cosh 2 (u>/2T) J 2^ [gcos(6» - a/2) - (1 + u/ji)].[qcas(0 + a/2) - (1 + cj/fi)] ^ 

where as above q = qkF/2mp, and the integration contour runs again infinitesimally above the real uj axis. Here in the 
angular integration we have taken the reference axis bisecting the angle between q! and q 2 . This angular integration 
is performed by residues, taking exp(i#) as a variable. For Y = (1 + u/jl)/q this leads to: 



" d A I = I (23 ) 

o 27r[cos(0-a/2)-Y].[cos(0 + a/2)-Y] (Y 2 - cos 2 (a/2))V^ 2 - 1 

where the cut in the determination of the square root has to be taken on the positive real axis, as it is clear when 
one considers the case of very large \Y\. When this result is inserted in Eq.(22) and Y is taken as new variable, one 
finds, introducing again the reduced temperature t — T/fi: 

/OO 1 y 
dY 5 o ; (24) 
.oo cosh 2 ((gy-i)/2i)(r 2 -r 2 )VF^T v 1 

with Yi = cos(a/2). 

Up to now we have made no approximation in our calculation and the result is valid at any temperature. Let us 
now focus on the low temperature regime t « 1 Because of the factor cosh~ 2 ((gY" — l)/2t), only the vicinity of 
Y = 1 will contribute. Moreover only the half circle contour around the pole Y = Y\ and the range Y > 1 contribute 
to the real part. In this last domain we have (qY — l)/2t > (q — l)/2t 3> 1 from our result for q Eq.(ll), so we 
can again simplify the hyperbolic cosine into an exponential (although this is not in practice a good approximation 
numerically). With the further change of variable Y = 1 + tv 2 /q in the resulting integral, we find to leading order: 

16tfl 2 J(a) = , =- — ; — — / dv P ' (25) 

P V ' acosh 2 [(l/4)ln§-/3 2 /2] V5F(1 + cos(a/2)) J v 2 + (3 2 V 1 

where we have substituted explicitely the result Eq.(ll) for q and have set 1 — (1 — cos(a/2))/i. The integral in 
this result can not be further simplified in general and is related to parabolic cylinder functions. Let us consider now 
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some important limiting cases for this result. First we take at fixed a the limit T — > 0. This implies /3 2 — > oo, the 
first term goes to zero and the integral is easily calculated in this limit, leading to: 

J («) = —TT ■ 2] ln s (26) 
4/x 2 sin 2 (a/2) 

On the other hand if at fixed T we take the limit a — > 0, we have /3 2 — > 0. The limiting behaviour of the integral is 
7r/(2/3) — 7T 1 / 2 as can be obtained through a by parts integration, the dominant divergent contribution from the two 
terms cancels out and we are left with [19]: 

'M - £i <*) 

which goes naturally to infinity for T — > 0. These two limits can be obtained more rapidly by making the proper 
simplifications from the start of the calculation Eq.(22). 

The two limiting cases which we have just considered show that, at low temperature, J (a) has a quite remarkable 
behaviour. For most of the range it is negative as it can be seen from Eq.(26) and it even goes to large negative 
values when a gets very small. This tends to favor states with small angle between wavevectors, as we will see below. 
On the other hand for a = or very small J (a) is positive and very large, as it results from Eq.(27). Clearly at low 
T the interesting range is the small a domain where we can write (3 — a/ (St) 1 / 2 . Surprisingly J (a) first starts to 
increase strongly from its a = value before going down to very negative values. This can be seen simply by looking 
at the specific point 1 = (1/2) ln(7r/2£) for which the second term is negligible and which gives to dominant order 
J (a) = 7r/(32^ 2 )t- 3 / 2 (ln( 7 r/2i))- 1 / 2 , even more diverging for T -> than the T = value. The integral of J (a) 
over a can also be analytically evaluated, and it shows that the strong positive peak at small a dominates over the 
negative contribution from the rest of the range. This quasi-singular behaviour of J (a) is summarized in Fig. 2 where 
we have plotted J (a) for various temperatures. 

We perform now the same kind of treatment for J(a). One goes again from Eq.(21) to an integration on the real 
frequency axis. However there is no integration by parts, and the angular average to be calculated is somewhat more 
complicated. It can nevertheless be performed by the same method and gives cxplicitely, with the same variable 

Y = (1 +(jj/p)/q as above: 

2 " <m_ 1 1 1 1 

o 27r[coa 2 (0-a/2)-Y 2 }.[cos{0 + a/2) 2 -Y 2 }^ 2YVY t ^T { Y 2 -cos 2 (a/2) + Y 2 -sin 2 (a/2y ( ' 

The second term is obtained from the first one by the change a — > n — a. This leads to J {a) = J (a) + J (ir — a) 
with: 

f°° 1 1 

8fl 2 q 2 J (a) = -Re J JY tanh((^ - l)/2t)-j== y2 cog2(a/2) (29) 

Here we have contributions coming from half circles around the poles at Y = ±cos(a/2) = ±Y\ and contributions 
from the two domains Y > 1 and Y < — 1. Neglecting terms which are exponentially small in the low temperature 
limit, we can replace the hyperbolic tangent by —1 in the Y < —1 domain and for the Y = —Y\ pole. Gathering 
similar contributions this gives: 

n2 2r,N t , „ ,..„, f°° „, 1 -tanh((gF- l)/2i) 2a - tt 

8p?q 2 J (a) = tanh gYi - l)/2i) + / dY ^ ' ' + — (30) 

V ; sina uy Ji (Y 2 - cos 2 (a/2))VY 2 ^T sina V ' 

where we have used the fact that, without the tanh((<jY" — l)/2i) term, the integral can be performed exactly to give 
(ir — a)/ sin a. We have taken advantage of this to have the factor 1 — tanh((gF — l)/2t) which is going rapidly to zero 
for large Y. This will be of use when we consider below temperature corrections. The last term in Eq.(30) disappears 
in the combination J (a) + J (tt — a), so we omit it from now on. 

Looking now for the dominant contribution at low temperature, we simplify the hyperbolic tangent in the range 

Y > 1 since its argument (qY — l)/2t is large and positive. So one finds expressions which are similar to the one 
encountered in the calculation of J (a). Moreover, because of the relation J(ir — a) = J (a), we can restrict ourselves 
to the case a < tt/2. In this case the expression obtained by the replacement a — > tt — a is very simple since the 
corresponding value of [3 2 = (1 — cos(a/2))/i is always very large at low T. This leads finally to: 

o 2 t, n 2tt ,„ 1 2 f°° , exp(-w 2 ) 

8fl 2 J(a) = — : (1 . .-.„ — -) + -,= dv p ' 31) 

r y J sina l + (2i/ 7 r) 1 /2 e xp(/3 2 ) ; ^ J v 2 + /3 2 v ; 
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where we have taken into account that the second term is only significant when a is small, so we have made a = in 
its prefactor. 

In the limit T — > with fixed a, we have (3 2 — > oo and one gets merely: 



K«) = - — — (32) 



1 

4/x 2 sin a 

Again if we take a — > at fixed T, /3 2 — > 0, the dominant divergent contribution cancels and we find: 



J (a) = -± (33) 

These two limits can be obtained again more directly. From these cases we can guess that J (a) is always negative. 
This is seen on Fig. 3 where J (a) has been plotted. However, in the same way as J (a), it has also a singular behaviour 
at small a. While for (3 ~ a/(8t) 1 ^ 2 » 1 it diverges as — 7r/(4^ 2 a), it goes to the finite value — l/(4^i 2 ) for a = 0. We 
note that the divergent behaviour in or 1 is weaker than the one in a~ 2 found for J (a) . Similarly J(0) » \J(0)\ at 
low T. So J(a) will play the dominant role and J (a) will only give a subdominant contribution. 



B. Temperature corrections 

The expressions Eq.(24) and Eq.(31) found above respectively for J (a) and J(a) are the leading terms at low 
temperature. While containing the dominant physical behaviour, we do not expect them to be so good quantitatively 
at intermediate temperature. It is actually possible to improve these analytical results markedly in this respect at the 
price of a slight complication by including the first two terms in an expansion in powers of t 1 / 2 . This is the similar to 
what we have done at the end of section II for the second order term, and we follow here the same procedure, together 
with the steps we have just taken above for the calculation of J (a) and J (a). 

We first consider first J (a) and start from the exact Eq.(24). Since the contribution from the domain Y < — 1 
is again exponentially small in the low temperature regime, we keep only the half circle contour around Y\ and the 
Y > 1 domain. The integral appearing in Eq.(24) can be split in two terms by: 



Y 1/1 1 



Y 2 - Y 2 2\Y-Y! Y + Yi 



(34) 



Therefore J (a) can be written as J {a) = J (a) + J (2ir — a) and we concentrate on the calculation of J . With the 
same notations a = (q — l)/2t, i — t/q and the change of variable (qY — l)/2i = a + u 2 , Eq.(24) leads to: 



32tfi 2 q 2 J (a) = -t- 1 / 2 i2 ^ 9 ~— 7^,1 r^~^ + , , ^ .,,"2, ^7 ( 3 5) 



di 

cosh 2 (u 2 + a) {u 2 + /3 2 /2)Vl + tu 2 ' sin(a/2) cosh 2 (a - /3 2 /2) 



where we have set j3 2 = (1 — cos(a/2))/t. The first term in the r.h.s. of Eq.(35) is expanded for low t by (1 + tw 2 ) -1 / 2 ~ 
1 — tu 2 /2 and the resulting temperature correction is calculated to lowest order, replacing a by a — ln(7r/2i)/4. This 
gives our following final expression, which has to be used together with Eq.(16) for a and q: 

«„- 2 2 T / n - 1/2 f°° du 1 TT 21 f°° , v 2 cxp(-v 2 ) 

mri 2 q 2 J (a) = -r 1 ' 2 / = =— - + = = + -= / dv , v - ; (36) 

^ H oy 1 J cosh 2 (u 2 + a)u 2 +(3 2 /2 sin(a/2) cosh 2 (a - /3 2 /2) Jo « 2 + P 2 

The second term in the r.h.s of this equation is only relevant for angles that are close to zero. Otherwise, it gives an 
exponentially small contribution in the low temperature limit. In the expression obtained by the change a — > 2n — a, 
this term is exponentially small in any case and can be forgotten. The first and the third term can be calculated 
more explicitly for angles a which are not close to zero (this is in particular the case for J (2ir — a)), which implies 
j3 — » +00. For the first term in particular we can make use of our results in section II for the temperature expansion 
where similar terms were found. This leads finally in this regime to: 

,„99 T /x 2 _2 — cos Oil 2 , . 

lSy 2 ? Jo (a) = — -— + t r-fc-f 37 

1 — cos(a/2) (1 - cos(a/2)) 2 

For various temperatures, we compare in Fig. 2 these low temperature expressions with the exact calculation of 
J (a), from the direct numerical summation in Eq.(20) over Matsubara frequencies (with q obtained by the numerical 
minimisation of I(q,j2,T) given by Eq.(6)). We see that, at this level of accuracy, they agree remarkably well up to 
rather high temperatures. 
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FIG. 2. J(a) for various values of the reduced temperature i = T/(qjl) = 0.01, 0.05 and 0.1. The lines are calculated using 
our low temperature expressions for J(a) and a, Eqs.(36) and (15), while the points correspond to the exact result numerical 
summation over Matsubara frequencies. The dashed line is the asymptotic behaviour for J (a) in the T — ► limit Eq.(26). The 
dashed-dotted lines are the leading order low temperature result Eq.(25) 



We turn now to a similar calculation for J (a) and take up from Eq.(30). We have already pointed out that, since 
we can restrict ourselves to < a < tt/2, the replacement a — ► 7r — a gives a value of /3 2 = (1 — cos(a/2))/t which is 
always large at low temperature. Therefore, in this replacement, the integral in Eq.(30) can be calculated to lowest 
order in t and gives 2tj cos 2 (a/2). In the same way in the first term the hyperbolic tangent can be replaced by -1 in 
this substitution. With the same change of variable as for J (a), we have finally: 



8fi 2 q 2 J (a) 



sin a 



(1 - tanh(a - (3 2 /2)) 



2t 



cos 2 (a/2) 



r c 

-1/2 / 

Jo 



du 



1 — tanh(u + a) 



(1 + cos(a/2) + 2tu 2 ){u 2 + /3 2 /2)Vl + tu 2 

(38) 



As for the calculation of J (a), we can have a more explicit result by expanding the denominator in the integral to 
first order in powers of i and computing the resulting temperature correction to lowest order by replacing o by s - 
But we will not write the cumbersome resulting formula. 

Nevertheless, as done previously for J (a), we compare in Fig. 3 our low temperature expressions with the exact 
calculation of J(a), q and p. The discrepancy appears only above i ~ 0.05 in the vinicity of a = tt/2 and is mainly 
due to the fact that it becomes inaccurate to consider that fi 2 is large in calculating J (tt — a), as we have done (see 
above Eq.(31)). However our low temperature expressions are already sufficient to describe the cascade. 
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FIG. 3. J{a) plotted for various values of the reduced temperature i — T/(qp.) = 0.01, 0.05 and 0.1. The lines correspond 
to the low temperature expressions for J(a) and a, Eqs. (38) and (15). while the dots show the exact results from numerical 
summation over Matsubara frequencies. The dashed line is the asymptotic behaviour Eq.(32). 



IV. CASCADE OF ORDER PARAMETER STRUCTURES 
A. Ingredients responsible for the cascade 

The second order term in the expansion of the free energy fl gives the location of the transition jl(T) and the 
optimum wavevectors modulus q entering the order parameter A(r) but it can not distinguish between different order 
parameter structures, because from Eq.(2) the Fourier components A g of the order parameter A(r) are decoupled in 
the second order term. The resulting degeneracy is lifted by the fourth order term when one goes slightly into the 
supcrfluid phase. This is the basis of the analysis of Larkin and Ovchinnikov [2] and of Shimahara [15]. So to speak, 
the wavevectors directions are independent at the level of the second order term while the fourth order one provides 
an effective interaction between these directions. Since the expression Eq.(17) for the fourth order term depends only 
on the angle a„ jm between the wavevectors q„ and q m through J(a n ^ m ) and </(a„ jrra ), the wavevector interaction 
appears therefore simply as pair interactions depending only on the relative positions of the 2D wavevectors on the 
circle. We note however that J(a n ^ m ) corresponds to an interaction between two pairs of opposite wavevectors (q, — q) 
whereas J(a n ,m) gives an interaction between two single wavevectors. 

We consider now the ingredients which lead to the prediction of a cascade of transitions between order parameters 
with increasing number of wavevectors when the temperature goes to zero. First, since J(a) takes very large positive 
values for a smaller than a critical angle a , we can consider that the angle domain < a < a is forbidden in order 
to avoid a dramatic increase of the free energy Q. To be quite specific we define ao by J(a ) = 0. However we could 
as well take the critical angle as a c which gives the minimum of J (a), since ao and a c are anyway very close as it 
can be seen on Fig. 2 (we will indeed consider also a c in the following). On the other hand it is favorable to take a 
slightly above ceo since it is the region where J(a) takes its most negative values. We note that the behavior of J(a) 
is not as strong compared to J(a), so we neglect J(a) in a first approximation. Next we see that J(0) diverges as 1/T 
at low temperature and it is necessarily present in the free energy Eq.(17) from the terms with q„ = q m . However 
their unfavorable effect on ft is lowered in relative value if one increases the number ./V of wavevectors, since we have 
N terms in Eq.(17) containing J(0) compared to a total number of N 2 terms. Since we want to minimize tt, this 
leads us to increase N and hence to decrease the angle between wavevectors as much as possible. Since the angle 
between wavevectors is bounded from below by a 0; we expect for symmetry reasons that the optimum structure to 
have regularly spaced wavevectors with an interval angle slightly above ao- The last ingredient to predict the cascade 
is the fact that ao decreases with temperature. As a result, the number N of plane waves in the optimum structure 
increases when the temperature goes to zero. This leads to a cascade of states where the T = limit is singular. In 
the next sections, we present our arguments in more details. 



B. Study of the cascade 



We begin by considering only the contribution Q a to the free energy ft from the J(a) terms. Even so the full 
problem of finding the order parameter structure minimizing the fourth order term is not a simple one. However it is 
natural to assume that the wavevectors of the N plane waves have a regular angular separation with an angle 2n/N 
between neighbouring wavevectors, so that their angular position is given on the circle by a n — 2m: /N. Otherwise 
we would have a minimum corresponding to a disordered situation for the angles, which sounds quite unlikely (note 
that we can not collapse an angle to zero since we work at fixed N; we will later on minimize with respect to N). 
Then it is easy to show that the weight w n = | A qn | 2 of the various wavevectors are all equal. Indeed the fourth order 
term Eq.(17) is just a quadratic form in w n and minimizing it is formally identical, for example, to find the lowest 
energy for a single particle in a tight binding Hamiltonian on a ring, with hopping matrix element J{a n ^ m ) between 
site n and site m (except for the on-site term which is J(0)/2). The eigenvectors are plane waves and the eigenvalues 
are J(0)/2 + ^2 n Z\ J{ot n ) cos(n</>fc) with cf>k — 2kn/N and k = 0, 1, ..,N — 1. Since J(a n ) < for n ^ 0, the lowest 
eigenvalue corresponds to k — which means that the weight w n are all equal. 

Conversely if we assume from the start that the all weights w n are equal, our problem of finding the best a„'s 
is the same as the one of finding the equilibrium position of atoms on a ring with repulsive short range interaction 
(because J(a) is large and positive for small a) and attractive long range interaction (because J (a) is negative for 
larger a). We expect the equilibrium to correspond to a crystalline structure with regularly spaced atoms. This takes 
into account that J (a) is a long range potential (clearly this regular spacing would not be the equilibrium if we had a 
strongly short range potential enforcing a specific distance between our atoms). Naturally in this argument we take 
N large enough to fill up the ring with atoms. So we come to the conclusion that the minimum energy corresponds 
to equally spaced wavevectors. 

Finally when we minimize the total free energy Eq.(2) with respect to the weight w n we find the general result: 

n_ » 2 (g,M,r) 2 

N ~ G 2 (N) lJyj 

where we have, in the case = Q a : 

N-l 

NG 2 (N) =2J(0) + 4^ J(a n ) (40) 

n=l 

We note that we still have a degeneracy of the lowest energy configuration with respect to the choice of the plane 
waves phases. 

We now take also into account in the free energy the terms containing J(a„ jm ). These terms appear when there 
are pairs of opposite wavevectors (q, — q). With our assumption of regularly spaced wavevectors, this corresponds to 
even N - in which case we have N/2 pairs - whereas for odd N, there are no J{a n ,m) terms in the free energy. Note 
that J(a) is negative for any angle a so that it is favorable to take pairs of opposite wavevectors. This seems to be 
in favor of taking N even and we will indeed see that as a result states with even N will be selected. 

For even N, the total free energy is: 

Q O N/2-1 N/2-1 

W = W + 2 ^2 E (l-^, m )A„A_„A^Al ro JK_ m ) (41) 

n=0 m=0 

where we have used J(tt — a) = J(a) and a n ^ m = 2(n — m)ir/N = a n - m ; A„ is a shorthand for A qn and so on. The 
fact that J (a) is always negative has a direct consequence on the phases of the plane waves. In order to minimize the 
free energy they have to be chosen so that A„A_„A^ l AI m is always real. Writing A„ = |A„|e^™, this implies that 
4> n + (f>- n = <&() for any n, where $o is a constant phase which can be chosen to be zero, since this merely corresponds 
to a global phase change for the order parameter. This link between phases for opposite wavevectors removes only a 
part of the degeneracy, since the phase <f>„ = —<t>- n can be still arbitrarily chosen. Now we see that the contribution 
of the J(a n ,m) terms to the free energy is a quadratic form in w' n — A„A_„ with negative coefficients, and no on-site 
contribution. This is quite similar to what we have found for the J (a) terms in the preceding subsection. So the 
equal weight distribution, which minimizes the J (a) terms, gives also independently the minimum of the contribution 
of the J(a) terms. This means that the inclusion of these last terms in the free energy does not change the optimum 
structure apart from the phase link between opposite wavevectors. Finally the total minimum free energy is still given 
by Eq.(39), where Eq.(40) is valid for odd N, but has to be replaced for even N by: 
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N-l N/2-1 

NG 2 {N)=2J(0)+A^J(a n )+A £ J(a n ) 

n—1 n—1 



(42) 



In this last case the corresponding equilibrium order parameter is real, being a sum of cosines of the form A(r) = 
| Ai| J2i cos(qi.r + <£>i). We turn now to the minimization of the free energy with respect to N. 




0.05 0.1 0.15 0.2 

i 

FIG. 4. [G2] 1 as a function of the reduced temperature i = T/qjl, for different values of N. Full lines correspond to our 
analytical low temperature expressions from Eqs. (15), (36) and (38), while dashed lines show the exact result. They agree 
perfectly for N = 4, 5, 6 in the i < 0.05 domain. 



C. Minimization with respect to N. 

We consider now the numerical calculation of G2(N) for different order parameter structures, i.e., for various 
values of the number of plane waves N. From Eq.(39) the equilibrium order parameter corresponds to the maximum 
[G2(N)] _1 . In Fig. 4 we show for [G2(N)}~ 1 both our low temperature expansion and the exact evaluation of the 
Matsubara sum in Eqs. (20) and (21), as a function of the reduced temperature t. We see that below t ~ 0.05, the 
low temperature analytical expressions, Eqs. (15), (36) and (38), agree remarkably well with the exact result. They 
are therefore completely sufficient quantitatively to study the cascade. 

Fig. 5 displays the same quantity for lower temperatures, exhibiting the cascade of transitions. An interesting 
feature of this cascade is that an order parameter with an odd number of plane waves is never the lowest energy 
solution, except from the N = 3 case, already found by Shimahara [15] . The reasons for this result will appear clearly 
in the next section, where we will study analytically the asymptotic regime of low temperatures. 

V. LOW TEMPERATURE ASYMPTOTIC BEHAVIOUR 
A. Minimum angle 

In the low temperature limit t <C 1, we can derive an explicit expansion for the value of the zero ao of J {a), as well 
as the location a c of its minimum. These critical angles are crucial to study the cascade since they give essentially 
the minimum angle between two wavevectors. In this low temperature range we have q ~ 1, so that t ~ t. As can 
be seen in Fig. 2, J(a) has two extrema: a maximum for a around a ~ /3 2 /2, as it is clear from the cosh -2 term in 
Eq.(25), and, for a slightly larger value of a, a minimum at a c . The zero ao is naturally in between. The condition 
/3 2 /2 ~ a ~ ao = (1/4) ln(7r/2i) implies a ~ ln(7r/2t). Therefore at low t, both ao and a c are in a domain where 
a — * and x = /3 2 /2 = a 2 /16< — > +00. In this regime the two terms in Eq.(25) for J(a) simplify to give: 
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FIG. 5. G2 plotted for temperatures lower and values of TV higher than in Fig. 4, using Eqs. (15), (36) and (38). For clarity 
we do not represent solutions of N larger than 14 but they can be easily numerically evaluated. 



fi 2 J{a) 



1 



I6ta cosh 2 X a 2 



(43) 



where X = x — clq. The second term in this expression simply comes from — l/4sin 2 (a/2) which is the asymptotic 
limit Eq.(26) for p, 2 J(a). 

Let us first consider the calculation of a c . The first term is responsible for the strong upturn near the minimum. 
From the derivative of Eq.(43) J (a) is extremal for: 



cosh 2 ^-* 



3/2 



2tl/2 (tanhX + -) 



(44) 



The minimum we are looking for is in the domain X > and from the above equation it is found for a large value of 
X. This allows to simplify this equation into: 



kx 3 ^ 



(45) 



with k = (2-7T 3 /i 2 ) 1 / 4 . Writing this equation x = In k+ (3/4) \nx one can generate a solution by the recurrence relation 
x n+ \ = lnfc + (3/4) lnx„ which converges very rapidly. For example the second iteration gives: 



x = lnk+ (3/4) ln[ln k + (3/4) ln(ln fc)] 
The corresponding result for a c is to leading order in the low temperature limit: 

(2tt 3 ) 1 /2^ 1 



8t In- 



t 



(46) 



(47) 



but numerically this is not such a good result at low temperature and one has rather to perform a few iterations to 
get the correct answer. For example in Eq.(45), the exact result is x = 5.94 for k — 100, while Ink = 4.60, but the 
second iteration Eq. (46) gives x = 5.92. 

Then in order to find the optimum number of plane waves we notice that J (a) rises very rapidly below a c , so we 
can not have basically the angular separation 2tt/N between two plane waves less than a c . Since on the other hand 
it is energetically favorable to take N as large as possible, as long as J(a) is negative, we can find an asymptotic 
estimate of the optimum value iV of iV by taking the integer value of 27r/a c , that is: 



N n = E 



2tt 

a,* 



(48) 



Following the same procedure, we can also derive an explicit expression for the angle «o corresponding to the zero 
of J (a). From Eq.(43), we find: 
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(49) 



Naturally a is close to a c since J(a) is rapidly increasing below a c . We are still in the X > domain and the 
solution corresponds to X large. This simplifies the above equation into: 



V2 



(50) 



where k has been defined previously. The corresponding recurrence relation which generates the exact solution is 
x n+ i = ln(fc/\/2) + (1/4) lna;„. The second iteration gives here: 



, k 1 i 
x = In — = ^ — In 

V2 4 



In —= + - In In —= 

V2 4 V \/2 



(51) 



and ao is simply given by ao = V / 16te. We can then make the same argument as above for a c , and write the following 
asymptotic estimate of the optimum value iV Q : 



N n = E 



2tt 
a 



(52) 



which will naturally be found to be quite near the above one Eq.(48), since ao and a c are quite close. 

Naturally we can calculate numerically exactly the optimal value -/V of N as a function of the temperature, from 
our above results. In the same process we find also the critical temperatures where the No changes. They are given by 
the crossings of the G*2(iV) curves for different values of N, as it is seen on Fig. 4 and 5. In Fig. 6 these exact critical 
temperatures are compared with the exact calculations of 2ir/a c and 27r/a 0j as well as their asymptotic values. As it 
is seen on this figure, it happens that the optimal values of N falls essentially just between 2it/a c and 2ir/ao in the 
low temperature domain. 




t 

FIG. 6. Optimal value of iV and locations of the critical temperatures. The dots give the temperatures where the optimal N 
changes: the white dot is the optimal N above the critical temperature whereas the black dot is the optimal N below. 2n/a c 
(curves very near the white dots) and 27r/o?o (curves very near the black dots) are also shown on this figure. The dashed lines 
show the exact calculations,while the solid lines are determined with the recurrence formulas below Eq. (45) and Eq. (50) with 
only one iteration. In the insert for clarity we only show the "devil" staircase corresponding to the white and black dots. 



B. Asymptotic evaluation of G2(N) 



It can be seen on Fig. 2 that J(a) switches rapidly above a c to its large angle asymptotic behaviour: 

1 



4 sin 2 (a/2) 



(53) 
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so we may with a fair precision use this simplified expression for J(a n ), with a n = 2im/N , to evaluate the sum in 
Eq.(40) for G 2 {N). At low temperature the number of plane waves N is large and the sum is dominated by the terms 
with small n corresponding to small angles, for which we have sin(a/2) ~ a/2. Hence we have for the sum coming in 
Eq. (40) for G 2 {N): 

4^ a „2*n, 8 1 N /1\ 

n—l n=l v 7 

For even N, we have also to evaluate the contribution coming from the terms in Eq. (42) containing J (a). Again 
we can use for most angles the asymptotic expression Eq.(32) for J(2irn/N). Once more the dominant contribution 
comes from the small n terms, and the logarithmic leading behaviour is accurately given by the asymptotic expression 
of J (a) as: 

JV/2-1 N/4 
n—l n—l ' 

[VERIFIER lc n=l]Wc sec here cxplicitely that J (a) has a dominant role compared to J (a). Taking into account the 
J(0) contributions, the dominant behaviour of G2(N) can be summarized in the limit of large plane wave number N 
and low temperature as: 

M 2 G 2 (A) ~ ^ - j Afodd (56) 

fi 2 G 2 (N) ~ - y -\nN + cte Neven (57) 

These expressions correspond actually in Fig. 5 only to the rising part of {G 2 (N)]~ 1 on the low temperature side. 
The downturn of [G 2 (N)]^ 1 for higher temperature is due to contributions from the positive part of J (a) which are 
beyond our asymptotic approximation Eq. (53). Nevertheless we have also found above the critical temperatures for 
switching from a given value of A~ to the next one. From Fig. 5 it is seen that it is enough to plug these critical 
temperatures into the evaluation of the rising part of [G 2 (N)]~ 1 to obtain the free energy at the transitions. When 
we substitute accordingly in these expressions Eq. (56) or (57) the value N = 2n/a c for the optimum plane wave 
number we have: 

> 2G >w=^( i -W +o{Hvt)) (58) 

Since St/a 2 <~ l/ln(l/t) from Eq.(47), we find that G 2 (N) is always positive in the low temperature range (more 
precisely we have 1 — (8ir 2 t)/(3al) = 1 — n 2 /6x which is positive as soon as x > 7r 2 /6, i.e. from Eq.(45) exactly for 
t < 0.62). This means that the transition stays always second order (a negative sign would have implied a first order 
transition and a breakdown of our fourth order expansion). 

Finally it seems from Eq. (56-57) that it is favourable to have an even number of plane waves in order to take 
advantage of the additional J contribution. This can be confirmed by a more careful comparison between two 
consecutives values of N. Let us assume that is odd with N = 2n/a c and compare G*2(A r ) to G 2 (N — 1). By going 
from N to N — 1 we gain naturally the J term in Eq.(57), but we increase the first term due to the unfavorable of 
the J(0) term. Nevertheless: 

fl 2 [G 2 (N 1) - G 2 (N)} = 2N{ ^_ 1)t ln(A - 1) + 0(1). (59) 

where the dominant term becomes exactly with the help of Eq.(45): 

2 5 1 32 

fl 2 [G 2 (N - 1) - G 2 (N)] = -x(l - — ) + In a + - In — (60) 

which is always negative since its maximum, reached for l/x m — (4/5)(l — 2/tt 2 ), is -0.107. Therefore the asymptotic 
evaluation of G 2 (N) shows indeed explicitely that the only order parameters which appears at low temperature are 
those with an even number of plane waves, in agreement with our exact numerical results. 
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C. Phenomenological interpretation 

It is interesting to compare our results to the pairing ring picture explored by Bowers and Rajagopal [20] (BR) in the 
three dimensional case. BR have looked at the free energy expansion at zero temperature in 3D, extending the work of 
Larkin and Ovchinnikov. In the case of multiple plane waves in the order parameter, they have pointed out a simple 
physical interpretation of their results. Their picture is based on the Fulde and Ferrell study [1] of a single plane wave 
order parameter where they showed that one may separate the wavevector space in two complementary domains: 
the pairing region and the pair breaking region. In the case of a vanishing order parameter and considering only 
wavevectors at the Fermi surface, the pairing region is merely an infinitely thin circle. It is given by the intersection 
of the up and down spin Fermi surfaces, after shifting one of them by q. The total opening angle of this circle is 
thercforer given by tpo = arccos(l/g). In 3D and at T = 0, q = 1.1997 which gives V'o = 67.1°. For a non-zero order 
parameter amplitude A, the circle broadens and becomes a ring whose width is given by A//2. Therefore, the pairing 
region for a single plane wave is a ring drawn on the Fermi surface, whose center is on the plane wave direction q (the 
restriction to the Fermi surface is justified by the fact that only wavevectors close to the Fermi surface are relevant 
for pairing). 

In the case of multiple plane waves in the order parameter, BR showed that the intersection of pairing circles from 
different plane waves is energetically largely unfavored, the worst case being when two circles are exactly in contact. 
The maximisation of the pairing regions leads naturally to increase as much as possible the plane waves number iV 
with the constraint that they are no intersecting circles. This leads to nine circles on the Fermi surface. BR rather 
concluded that the least energy state was the one with eight plane waves for 'regularity' reasons. In particular, this 
leads to a real order parameter. We see now that our cascade in 2D is a direct consequence of these two heuristic 
principles : the no-intersecting circles rule and the maximisation of plane waves. 

In 2D, the circles are replaced by pairs of points separated by the angle ipo, and the no-intersecting rule becomes 
now a no- entanglement rule between different pairs. But this is just equivalent to our above finding from Fig. 2 that 
there is an effective short range repulsion between different wavevectors q, with indeed the worst case corresponding 
to the contact since the maximum of J(a) is just below its zero. The contact condition a — p 1 /2 that we have 
found (corresponding to the maximum of J(a)) can be rewritten as cos(a/2) = l/q which implies a/2 = ipo/2. This 
is exactly the same as the contact condition of BR. Next we have found explicitely that it is indeed favorable to 
maximize the number of plane waves, taking the short range repulsion into account. Our cascade is due to the fact 
that, when the temperature goes to zero, ipo also drops to zero and the pairing region for a single plane wave shrinks 
to a single point, corresponding to the situation where the two up and down Fermi circles are just in contact. 

VI. CONCLUSION 

In this paper we have investigated the low temperature range for the FFLO transition in two dimensions and we 
have shown that the order parameter is no longer a simple cos(q.r) at the second order transition, in contrast with the 
situation found near the tricritical point. Instead the transition is toward more and more complex order parameters 
when the temperature goes to zero, which gives rise to a cascade with, in principle, an infinite number of transitions. 
At the transition these order parameters are in general a superposition of an even number of plane waves with equal 
weight and equal wavevector modulus, which corresponds to a real order parameter equal to a superposition of cosines. 
The directions of these wavevectors are found to be equally spaced angularly, with a spacing which goes to zero when 
the temperature goes to zero, which is the reason for the ever increasing number of plane waves in this limit. The 
singular behaviour in this limit T = is actually present in all physical quantities. It arises because, in order to 
obtain the lowest energy, the two Fermi circles corresponding to opposite spins come just in contact when one applies 
the shift corresponding to the wavevector q of the FFLO phase. This is this situation, with these two circles just 
touching each other, which gives rise to the singularities. Naturally this is linked to the fact that in 2D the Fermi 
surface is actually a line. Hence this singularity is a general feature of 2D physics and we may expect it to give rise 
to similar consequences in more realistic and more complex models describing actual physical systems. 

Naturally in this paper we have addressed this remarkable situation in the most simple physical frame. We have 
considered the simplest BCS model with isotropic Fermi surface, that is a circle for our 2D case. Naturally if anisotropy 
is introduced in the dispersion relation of the electrons, or in the effective electron-electron interaction, or both, the 
physics will be more complex. In the same way we expect Fermi liquid effects [9] to bring important quantitative 
modifications. Similarly if we consider a quasi 2D superconductor with a magnetic field not perfectly aligned with the 
planes, currents will be produced and orbital contributions to the free energy will arise. However in order to have the 
physics right in all these complex situations, it is quite clear that one has to obtain the correct limit in the simplest 
possible case that we have considered in the present paper. 
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VII. APPENDIX A 



In this section we rederive rapidly our result Eq.(ll) for the leading order of the optimal wavevector at low 
temperature, taking the same starting point as Bulaevskii [18]. 
We start from Eq.ll of Rcf. [18]: 

1 f Pv F 1 

HT C /T) = -Re / dn Ml/2 + z(2^H + Q)/4ttT)) - ^(1/2)] (61) 

K J-Pv F \/P 2 V F - ft 2 

We go to our notations by setting H = p,, P = q and introduce as above our reduced variables q = qvp/2p and 
t = T/p,. When we perform the angular integration as in Ref. [18] by setting f2 = Pvf cos 9 and introduce for the 
digamma function the integral representation (Eq.13 of Ref. [18]): 



l ' 2X ~ 2 L (^ + ^)(!xp(2^)-l) 



(62) 



we obtain Eq.14 of Ref. [18] (correcting some minor misprints). We can then expand the result to all orders in t and 
show that all the terms in the expansion are zero. However rather than displaying the steps of this calculation here, 
it is more convenient to immediately remark that the optimal q is obtained by writing that the derivative of Eq. 1 1 
with respect to q is zero. This provides an equivalent calculation (we work on the derivative instead of working on 
the function), which is somewhat easier and allows us also to make in the following the direct contact with our result. 
This condition on the derivative, which gives the equation for the optimal q, is: 

/■7T jn 

Im — V>'(l/2 + i(l + gcos(9)/27rf)) = (63) 
Jo n 

Taking the derivative of the above integral representation, we obtain a corresponding representation for ip' on which 
the angular integration is easily performed (the integrals can be found in Gradstein and Ryzhik [21]). We display the 
result by introducing the function f(x) = xj \f\x 2 — 1) and setting xq = 1/q , e = —iirt/q and x = x + e. The result 
is: 

Re [1 - f(x) + ef'(x) - £ f° . t, , (/(*) - /(*o + e(l + 2iy)))] = (64) 
2 J-oo sinh (ny) 

where the integral goes from — oo to oo because we have collected two terms into one. Now we can write for the first 
three terms the Taylor expansion: 

oo n+1 

1 - f(x) + ef(x) = 1 - /(so) + £ —. (x ) (65) 

^ n! n+1 

When we take the real part, only the odd order derivatives of / contribute because q is near 1 with q > 1, which 
makes x < 1 and Ref^ 2p ^(x ) = 0. This gives: 



oo 



2p+l o 

Re [1 - f{x) + ef'(x)\ = 1 + E ^2pjl2^TT f {Xo) (66) 

Similarly we can perform the expansion in the integral and expand (1 + 2iy) n by introducing the binomial coefficients 
C%. All the resulting integrals can be found in Gradstein and Ryzhik [21], and expressed in terms of Bernoulli numbers 
B 2m - This leads to: 

2 Re / - m h 2 (n ) [fix) f{X ° + 6(1 + 2lV))] = S (27TT)! /(2P+1)(X0) S C 2p+i4 m B 2m (67) 
Now a remarkable identity (found for example in Gradstein and Ryzhik [21]) for Bernoulli numbers states that: 

J2 C 2 2 ™ +1 ± m B 2m = 2p (68) 



rn — l 
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As a result, by gathering all the terms, all the coefficients of the powers of e are zero, and it is not possible to satisfy 
the condition that the above derivative is zero. The answer to this puzzle is that the contribution of the terms we 
have considered is not exactly zero, but exponentially small, which explains why we find it to be zero in a perturbative 
expansion. This is shown now in the following. 

We start again from the above condition Eq.(64) found for the optimum wavevector. We transform the integral by 
shifting the integration contour toward the upper complex plane by i/2 for the variable y. First we have to take care 
that integrant in the above integral has no singularity for y = 0, because proper cancellation between various terms. 
Hence we can also say that this integral is equal to its principal part. Next if we want a complete contour C, we have 
to add to this principal part the contribution of an infinitesimal semi-circle around y = with positive imaginary 
part. This contribution is easily found by residues to be equal to —ef'(x), so we have: 

Re [ef{x) - f ^ . (f{x) f(x + £ (1 + 2iy)))} = -^Re / , - [fix) f(x + e(l + 2iy))] (69) 

* J-oo sinh (Try) I Jc smh (Try) 

Now we set ny = z + m/2 and make use of isinh(z + in/2) — zcosh(z) which gives: 

Re [ef'(x) - | f° . (f(x) f{x + e(l + 2iy)))] = Uc f°° -^-(f(x) f{x + 2tz/q)) (70) 

* J-oo smh (Try) ^ J-oo cosh z 

where the contour runs infinitcsimally below the real z axis. The first term in the last integral is just equal to fix), 
so our equation becomes: 

* J-oo cosh {iry) 

It is easily checked that this equation is identical to our Eq.(10). Note that the fact that the contour runs infinitesimally 
above or below the real axis is unimportant since we take the real part, and the contribution along the cut due to the 
square root is purely imaginary. 

In conclusion we have obtained our basic equation for the second order term by taking the same starting point as 
Ref. [18]. The end of the argument to obtain Eq.(ll) is naturally the same as following our Eq.(10). In particular this 
argument shows that the integral in the above formula is proportional to exp [— (g — l)/i], so it can not be expanded 
in powers of t. 

* Laboratoire associe au Centre National de la Recherche Scientifique et aux Universites Paris 6 et Paris 7. 
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